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Abstract. We apply the immersed boundary (or IB) method to simulate deformation and detach¬ 
ment of a periodic array of wall-bounded biofilm colonies in response to a linear shear flow. The 
biofilm material is represented as a network of Hookean springs that are placed along the edges 
of a triangulation of the biofilm region. The interfacial shear stress, lift and drag forces acting 
on the biofilm colony are computed by using fluid stress jump method developed by Williams, 
Fauci and Gaver [Disc. Contin. Dyn. Sys. B ll(2):519-540,2009], with a modified version of their 
exclusion filter. Our detachment criterion is based on the novel concept of an averaged equivalent 
continuum stress tensor defined at each IB point in the biofilm which is then used to determine 
a corresponding von Mises yield stress; wherever this yield stress exceeds a given critical thresh¬ 
old the connections to that node are severed, thereby signalling the onset of a detachment event. 
In order to capture the deformation and detachment behaviour of a biofilm colony at different 
stages of growth, we consider a family of four biofilm shapes with varying aspect ratio. Our 
numerical simulations focus on the behaviour of weak biofilms (with relatively low yield stress 
threshold) and investigate features of the fluid-structure interaction such as locations of maxi¬ 
mum shear and increased drag. The most important conclusion of this work is that the com¬ 
monly employed detachment strategy in biofilm models based only on interfacial shear stress 
can lead to incorrect or inaccurate results when applied to the study of shear induced detachment 
of weak biofilms. Our detachment strategy based on equivalent continuum stresses provides a 
unified and consistent IB framework that handles both sloughing and erosion modes of biofilm 
detachment, and is consistent with strategies employed in many other continuum based biofilm 
models. 
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1 Introduction 

The subject of this work is the flow-induced deformation of a biofilm colony, which is a mesoscale 
collection of bacterial cells held together by an extracellular polymeric network (EPS) that is se¬ 
creted by the cells. The dimensions of a biofilm colony can be anywhere from tens to hundreds of 
microns, whereas the size of an individual bacterial cell making up the colony is on the order of 1-5 
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microns; our focus is on continuum models that treat the biofilm as a viscoelastic solid continuum 
rather than incorporating the dynamics of individual bacteria. The flow-induced deformations 
of the biofilm colony affect the fluid dynamic forces acting on it, and thereby also alter both the 
extent and the mode of detachment (i.e., sloughing or erosion) that may be experienced by the 
biofilm. We are particularly interested in understanding whether biofilm colonies gain any protec¬ 
tion against detachment when they are in close proximity to other colonies. To this end, our aim is 
to develop a robust numerical method for simulating the interaction between a biofilm colony and 
the surrounding fluid that is capable of capturing the different modes of biofilm detachment. 

Our approach is based on the immersed boundary (or IB) method, in which the biofilm con¬ 
tinuum is replaced by a network of Hookean springs. Although the IB method has already been 
used by several authors for studying biofilm deformation and detachment [1,26], our approach of 
computing an equivalent continuum stress at each IB node is markedly different and using it to 
initiate detachment provides us with a way of handling detachment in a manner that is consistent 
with other continuum mechanics based models such as [13]. 

1.1 Bacterial biofilms 

Bacterial biofilms are aggregations of microbes that grow on surfaces in an aqueous environment. 
They form when bacterial cells suspended in the fluid attach themselves to a surface and begin 
producing an extracellular polymeric substance in which the growing bacteria cells embed them¬ 
selves. Biofilms play contrasting roles depending on the scenarios in which they are encountered: 
in waste-water treatment [39] and environmental engineering they play a helpful role; whereas 
biofilms encountered on medical devices or food-processing equipment can be detrimental in the 
sense that they cause infections [9] or compromise food safety [48]. 

In these and other applications, biofilms experience vastly different physical conditions (tem¬ 
perature and pH), hydrodynamics (ranging from creeping to turbulent flow) and chemical envi¬ 
ronments (rich or sparse in nutrients, or saturated with antibiotics). Their adaptability to diverse 
conditions is believed to derive at least in part from the mechanical and chemical protection pro¬ 
vided by the gel-like EPS layer. It is important from an engineering standpoint to understand the 
mechanical properties of biofilms so that we can devise effective methods for not only enhancing 
biofilm growth and survival but also removing them from surfaces. Consequently, over the past 
10 years significant effort has been expended to develop novel rheological methods for measuring 
mechanical properties of biofilms [25]. Such experiments have typically arrived at different con¬ 
clusions on how to characterize biofilms, with some concluding that they behave as elastic [40] or 
viscoelastic solids [28], while others liken biofilms to viscoelastic fluids [55]. Furthermore, even 
within the same material class, measured material parameter values can vary over a fairly wide 
range. The general consensus is that biofilms behave as (visco-)elastic solids when the applied 
fluid shear stress is at or near the stresses at which the biofilm was grown, while at higher values 
of shear stress the biofilm can yield and behave as a viscoelastic fluid. In this paper, we restrict 
ourselves to the low shear stress case and model the biofilm mechanically as a viscoelastic solid 
embedded within a viscous fluid. 

1.2 Mathematical models of biofilm growth, deformation and detachment 

Mathematical models for simulating biofilm growth must take into account a wide range of mech¬ 
anistic and other dynamical processes, including growth and death of bacteria, attachment of cells 
to the substratum, transport of solutes (nutrients, metabolic products, antibiotics) within the sur- 
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rounding fluid and the biofilm itself, reaction kinetics, and removal of bacterial cells as clusters 
(sloughing) or as individual cells (erosion) when the biofilm EPS matrix weakens in response to 
fluid shear or chemical treatments. The earliest models developed in the 1980's [29,61] assumed 
that the biofilm is one-dimensional, while many subsequent modelling efforts have attempted to 
capture the 2D or 3D morphology that develops during biofilm growth processes. The treatment of 
the hydrodynamics and its interaction with the biofilm has received varying degrees of treatment 
in these models, ranging from some models that consider biofilm transport processes in isolation 
and use a specified nutrient concentration boundary layer thickness at the biofilm-fluid interface, 
whereas other models perform full fluid flow and nutrient transport calculations with or without 
incorporating fluid-structure interaction effects. A variety of approaches have been developed to 
model the growth and spreading of the biomass, including individual-based models [32], cellular 
automata [64], continuum models [2,13,30] and phase field models [67]. An exhaustive review 
of the different modeling approaches can be found in the article [60] and performance benchmark 
comparisons are available in [14]. 

Among the more complete models are those that capture multi-dimensional growth and fluid 
flow [13, 15, 16,44], although the role played by hydrodynamics in inducing biofilm deformation 
and detachment has most often received only ad hoc or approximate treatment. These approxima¬ 
tions are aimed at capturing the solid mechanics governing the dynamics of the deforming biofilm 
while simplifying as much as possible the complex fluid flow that surrounds it. The first attempt 
at solving the solid mechanics problem inside the biofilm and using it to initiate detachment was 
in [44] where the biofilm was treated as a linearly elastic material and a von Mises yield stress 
criterion was used to initiate detachment; however, this work neglected the effects of biofilm de¬ 
formation. In contrast, the particle-based biofilm model in [65] ignored the fluid and proposed 
a method for initiating detachment using a detachment speed function that is based on the normal 
velocity of the biofilm-fluid interface. A similar approach has been used in [49], which includes the 
effects of both biofilm growth and flow by employing a detachment speed that depends on the in¬ 
terfacial shear stress. This approach to initiating detachment rests on the assumption that regions 
where the interfacial shear stress is highest correspond to locations where the biofilm strain is also 
high. 

One of the aims of the current study is to verify the validity of this last assumption by consider¬ 
ing a full fluid-structure interaction simulation that is capable of determining stresses in both fluid 
and biofilm. We do not explicitly model biofilm growth but rather mimic the effects of growth 
by considering a family of biofilm colony shapes of different aspect ratio, where each colony size 
corresponds to a different instant of time during the growth of the same colony. In each case, 
we investigate how the deformation affects both the mode and the extent of biofilm detachment. 
This is in contrast with other approaches [13,44] where the solid mechanics are simulated but any 
deformations of the biofilm are neglected. 

1.3 Fluid-structure interaction in biofilms 

During the last several years, significant progress has been made in the study of fluid-structure 
interaction (FSI) in biofilms, driven by the increased availability of experimental data on biofilm 
mechanical properties and the motivation to understand the role they play in biofilm survival. 
Most FSI studies neglect biofilm growth by taking advantage of a natural separation of time scales, 
in that growth processes are very slow in relation to fluid motion and biofilm deformation. More 
recently, a phase field method has been applied successfully in 2D [34] and 3D [47] to simulate 
biofilm growth coupled with deformation arising from interaction with the surrounding flowing 


4 


fluid treating the biofilm continuum as a multiphase polymeric gel. With the exception of these two 
works, the most common approach in the biofilm FSI literature combines a Lagrangian discretiza¬ 
tion of the biofilm with an Arbitrary Lagrangian Eulerian (or ALE) formulation for the fluid. The 
first study of this kind appeared in [54] where they studied the deformation of a 2D hemispherical 
biofilm colony placed in a turbulent flow using the ANSYS® commercial software package. More 
recently [53], a similar ALE approach was used to study flow-induced oscillations of 2D biofilm 
streamers and their effect on mass transfer. A more realistic 3D biofilm model was studied in [7] 
that used sliced 3D confocal laser scanning microscopy data to construct the colony shapes, and 
employed a nonlinear hyper-elastic constitutive model that accommodates detachment based on 
a von Mises yield stress criterion. 

The aforementioned approaches have the advantage of being well-established in the literature 
and capable of easily incorporating experimental parameters that measure biofilm rheology. The 
primary disadvantage is their high computational cost and algorithmic complexity that result from 
needing to constantly re-mesh the fluid domain as the biofilm colony deforms. Moreover, this re¬ 
meshing cost increases enormously if detachment is incorporated in the model. 

Motivated by the desire to develop a simpler and more efficient computational approach, Alp- 
kvist and Klapper [1] proposed an alternate FSI strategy based on the immersed boundary (or IB) 
method in which the biofilm is discretized at a set of moving Lagrangian points. These IB points 
move relative to an underlying fixed Cartesian grid on which the fluid equations are solved. The 
elastic properties of the biofilm are captured by distributing forces onto the fluid that derive from 
a network of Flookean springs joining the IB points. An incompressible fluid pervades both fluid 
and biofilm regions and the biofilm inherits the density and viscosity of the surrounding fluid, 
so that the biofilm is actually a visco-elastic composite material consisting of elastic spring forces 
and fluid viscous forces. This approach has the clear advantage that no re-meshing of the fluid 
grid is required. To illustrate the versatility of their approach, Alpkvist and Klapper simulated the 
deformation of both 2D and 3D biofilm structures, while also incorporating a simple detachment 
criterion based on cutting individual springs when they are stretched beyond a critical length. We 
remark that other IB models for biofilms were developed prior to [1], namely the work of Dillon 
and collaborators [10,11]; however, these authors were concerned with slow flow and the dynam¬ 
ics of individual bacterial cells aggregating and settling on the substratum, and hence the results 
are relevant to different phenomena occurring on much smaller spatial scales and at much earlier 
stages of biofilm formation. 

A number of other IB approaches have since appeared, such as [57] who performed 3D simu¬ 
lations of biofilm deformation, comparing Hookean (elastic, spring-only) and Kelvin-Voigt (visco¬ 
elastic, spring plus dashpot) models for the biofilm material. Using a parametric study of spring 
stiffness and damping coefficient and comparisons with experiments, they established that realis¬ 
tic biofilm deformation behaviour can be obtained using the IB method. More recently, Hammond 
et al. [26] developed more detailed 2D and 3D IB models for fragmentation of a biofilm colony in 
which the location of actual bacterial cells from 3D images was used to determine coordinates of 
IB points. They employed a similar spring network and detachment strategy as in [1], but they 
allowed the density of the biofilm to differ from that of the fluid, and in subsequent work [27] also 
extended their approach to handle variable viscosity. 

Despite the increasing popularity of the IB method in biofilm FSI studies, two main challenges 
remain to be addressed before the method is capable of simulating realistic biofilm deformation 
and detachment. The first relates to connecting values of the IB spring parameters (elastic stiff¬ 
ness and damping coefficients) to actual biofilm material properties. Despite the parametric study 
in [57] that showed it was possible to determine suitable parameters for given biofilm colony shape 
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and IB spring network topology, there remains as yet no a priori method for determining IB param¬ 
eters for a given biofilm. 

The second challenge relates to initiating detachment in the IB framework in a way that is 
consistent with other more established continuum-mechanics-based biofilm studies mentioned in 
Section 1.2. Whereas IB methods have so far used strain in any given spring as a measure for 
initiating detachment, this is not a true measure of strain in the continuum mechanics context. 
Indeed, Hammond et al. [26] demonstrated that when using spring strain as a detachment criterion 
the resulting biofilm shape is sensitive to the critical strain parameter, and so it is unclear how to 
choose this parameter to match a given set of biofilm mechanical properties. 


1.4 Objectives and outline 

In this paper, we aim to address the challenges identified at the end of the preceding section by 
developing an immersed boundary approach for simulating biofilms that is capable of capturing 
realistic deformation and detachment behaviours. We begin with a 2D IB model inspired by that of 
Alpkvist and Klapper [1], and extend this work guided by two main objectives. Our first objective 
is to develop a novel approach for initiating biofilm detachment that is consistent with methods 
employed in biofilm models using continuum mechanics based calculations to enact detachment. 
In this way, we can retain the simpler spring network representation of the biofilm continuum and 
the advantages it offers, while also implementing a more physically realistic criterion for detach¬ 
ment. Our second major objective is to investigate the effect of flow-induced deformation on both 
the mechanical stability and mode of detachment experienced by weak biofilm structures having 
realistic shapes that resemble those grown under mass transfer limited conditions. This will al¬ 
low us to better understand such fundamental questions as how flow-induced deformation affects 
the forces acting on biofilms, and how spatial clustering of biofilm colonies can shield them from 
detachment by reducing the hydrodynamic shear forces. 

Our modelling approach incorporates the work done in several previous computational stud¬ 
ies of two-dimensional biofilms in [52,66] wherein we investigated the fluid shear-induced de¬ 
tachment forces (drag and lift) acting on rigid biofilm colonies that are both uniformly and non- 
uniformly spaced. In contrast to these earlier studies that were restricted to values of shear rate 
exceeding 10 cm/s, we focus in this paper on more flexible weak biofilm colonies immersed in a 
slower shear flow having shear rate less than 1 cm/s. 

The organization of the remainder of this paper is as follows. In Section 2, we define the prob¬ 
lem geometry and parameters, and describe the governing equations and corresponding numeri¬ 
cal scheme for our basic IB model framework. Section 3 contains the novel biofilm-related aspects 
of our IB model where we derive our approach for calculating fluid shear stress along the biofilm- 
colony fluid interface and from that the drag/lift forces acting on the biofilm colony. In Section 3.1 
we describe how we implement a modified version of the exclusion filter devised by Williams et 
al. in [63] for accurately approximating interfacial shear stresses on curved immersed boundaries, 
and then in Section 3.2 we describe how we adopt the concept of a continuum stress around each 
IB node in the biofilm material. Section 3.3 discusses how these quantities are used to obtain a 
more realistic detachment criterion for biofilm colonies in the IB context. Finally, in Section 4 we 
perform a number of numerical tests that validate our numerical approach and also demonstrate 
the advantage of our IB model for simulating realistic deformation and detachment events in fluid 
shear-induced biofilm dynamics. 
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2 Immersed boundary (IB) method 

The IB method is both a mathematical formulation and a numerical method for simulating the 
complex interaction between a deformable solid material and a surrounding fluid. The approach 
dates back to work of Peskin [42] who originally developed the approach to simulate blood flow 
interacting with heart muscle, and more recent theoretical and computational developments are 
summarized nicely in the review paper [43]. The IB approach has been applied extensively to the 
study FSI in bio-fluid mechanics, including such problems as amoeboid locomotion [8], platelet 
aggregation [19], and sperm motility [18]. 

The IB method is a mixed Eulerian-Lagrangian approach wherein the fluid equations are dis¬ 
cretized on a fixed, rectangular (Eulerian) mesh whereas the elastic structure is defined by a set of 
(Lagrangian) IB points that move relative to the underlying fluid mesh as the structure deforms. 
The effect of the immersed boundary on the fluid is represented using a singular source term in 
the fluid momentum equations, which is distributed onto the underlying fluid grid by writing it 
as the convolution of an elastic force density with a regularized delta function. 

2.1 Problem geometry 

A diagram of the problem domain is given in Fig. 1, which depicts shape of a biofilm colony en¬ 
countered during different stages of growth. The domain is delimited by two horizontal walls 
separated by a distance H, where the bottom wall (serving as the biofilm substrate) is held sta¬ 
tionary while the top wall moves horizontally with constant velocity U wa u. The domain is filled 
with a viscous, incompressible fluid so that in the absence of biofilm, the top wall will generate 
a flow that over time approaches a steady-state corresponding to linear (planar) shear with shear 
rat eG = U waU /H. 

Biofilm colonies having identical shape and size are placed along the bottom wall with a uni¬ 
form spacing of D/, between adjacent colonies. In order to simulate biofilms at different growth 
stages occuring under mass transfer limited conditions, we select four representative biofilm shapes 
with increasing aspect ratio that have fixed width but increasing height. Motivated by shapes at 
early stages of growth observed in experiments [12,31] and 3D simulations [13,15,30,67], we choose 
an idealized biofilm shape described by the equation 

(w^) +(l|) =1 with n 1 ,» 2 >0 and (2.1) 

which is half of a super-ellipse with width TV/, and height H[,. In the case when Yi\ = «2 = 2, Eq. (2.1) 
reduces to half of a regular ellipse, while further constraining Wj, = 2Hj, yields a semi-circle with 
radius Hj,. 

We believe that this choice is a reasonable parameterization of the typical finger-like shapes 
that dominate earlier stages of biofilm growth, in contrast with later stages that exhibit a classi¬ 
cal mushroom-shaped colony having an ellipsoidal head on top of a thin cylindrical stem. For the 
purpose of this study, we take values of ?ii = «2 = 2.5 and W/, = 50 pm, and select a sequence of 
three initial shapes with H), = 25, 50 and 75 pm respectively. In addition, we consider two other 
shapes: a semi-circle with diameter 40 pm (that is, ri \ = ri 2 = 2, TV/, = 40 pm and H/, = 20 pm); and 
an irregular mushroom-like structure with Wj,^60pm and FT/, ?=280 »m, whose shape is extracted 
from figures in [1], The latter mushroom shape is the tallest biofilm colony depicted in Fig. 1 and 
is chosen to illustrate the effectiveness of the detachment criteria developed in this study by direct 
comparisons with the results from [1,26]. 
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Figure 1: The problem domain consisting of a stationary bottom wall (e-f) separated by a distance H from a top wall (g-h) 
that moves at velocity (U wa u, 0). Biofilm colonies of width Wj, and increasing height represent shapes at different stages 
of growth. The walls are treated with IB points and the whole flow region is embedded in a slightly larger computational 
domain (abed) that is periodic in both x and y. 

Rather than imposing wall boundary conditions directly on the fluid, we simplify the fluid 
solver by treating the top and bottom walls using immersed boundaries and embed the problem 
domain inside a slightly larger fluid domain denoted by dashed lines in Fig. 1. Because we are 
interested in studying repeating arrays of biofilms colonies, periodic boundary conditions are im¬ 
posed in both the x and y directions so that the computational domain is O = [0, W] x [0 ,H] with 
W = Wb + Djj (labeled abed in the figure), and it contains a single biofilm colony. 

2.2 Governing equations 

The domain Q is filled with a Newtonian incompressible fluid having constant density p (g/cm 3 ) 
and dynamic viscosity y (g/ cms). Denote the infinitesimally thin immersed boundaries represent¬ 
ing the fixed bottom wall and moving top wall by Tand T tnf , respectively, and let Trepresent 
the solid elastic structure corresponding to the biofilm colony. Note that T^ 0 is actually a composite 
material that consists of the elastic force-generating material and fluid that co-exist within the same 
region. We also assume for the purposes of this study that the biofilm is neutrally buoyant and has 
the same density as the fluid (although it is straightforward to extend the IB approach to deal with 
variable density problems [22,26]). The fluid is therefore governed at all points x= (x,y) G O by 
the incompressible Navier-Stokes equations 


3u 7 

p-+pu- Vu=-Vy+yVu+f. 

(2.2) 

V • u = 0, 

(2.3) 
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where u(x,f) (cm/s) is the fluid velocity and p(x,f) is pressure (g/cms 2 ). 

The effect of the solid boundaries on the fluid is encompassed in the fluid forcing term f, which 
we consider next. Assume that the configuration of the solid material making up both channel 
walls and biofilm is described by a function X(q,t), where q is a generalized (dimensionless) pa¬ 
rameterization that is either a scalar (q = s) in the case of the walls T& of and T top , or else a vector 
(q = (r,s)) for a solid region like r/„ r) . The IB force in both cases is specified in terms of a discrete 
network of IB points connected by springs, and more detail on the precise form of these spring- 
force connections for walls and biofilm is provided later in Sections 2.4 and 2.5. Assume that the 
force generated by any deformed configuration can be described by an IB force density function 
F(X(q,f) depending on the current stretched configuration of the spring network. Then the fluid 
force f may be determined by spreading the force density at IB points onto the fluid using a delta 
function convolution 

f(x,f) = ^F(X,f) 5(x-X(q,t))dq, (2.4) 

where S(x) = 5(x)5(y) is the Cartesian product of two ID Dirac delta functions and T = TIJ U 
T/, of represents the set of all immersed boundaries. 

The final equation required to close the system is an evolution equation for the immersed 
boundaries, which we assume move with the same velocity as the surrounding fluid 

— = / u(x,f)^(x — X(q,t)dx. (2.5) 

at J n 

This is simply another way of stating the no-slip condition for a deformable boundary at location 
X(q,f). 

2.3 Numerical algorithm 

We now describe the basic numerical algorithm for solving Eqs. (2.2)-(2.5), which is a semi-implicit 
scheme very similar to the one employed in [50]. The fluid domain Q is discretized on a regular 
Cartesian mesh with coordinates x t j = (.r,,i/ ; ) = ( ih x ,jh y ) for i = 0,1,..., N x and j = 0,1,...,N ;/ , where 
h x = W / N x and hy = H/ Ny are constant grid spacings in the x and y directions (and we assume for 
simplicity that h x = hy). The time interval of interest [0,T] is likewise divided into a sequence of Nt 
equally-spaced points denoted t n = nAt for n = 0,l,...,N t , where At = T/N t is the time step. Let the 
discrete values of the velocity and pressure be denoted by u". and p'' respectively. Suppose that 
the immersed boundary is described by a set of IB points whose locations at any time t„ are 
given by X" = (X",Y") for l = 1,2.. .,N/,. The corresponding force densities are denoted by F”, with 
the precise specification of the immersed boundary discretization and force density calculation F" 
being given in the the following two sections. 

We now describe our algorithm for updating fluid grid quantities u "■ and p” and the IB config¬ 
uration X'J from time t n to time t n+ j. The algorithm proceeds in four main steps: 

Step 1: Compute the force density F" based on the current IB configuration X” as described in 
Section 2.5. 

Step 2: Spread the force density onto fluid grid points using a discrete representation of the delta- 
function convolution in Eq. (2.4) 

N b 

f//=E F ?4(Xi;-X?)A 

i =i 


( 2 . 6 ) 


9 


where Sh( x ) is a regularized delta function given by 


with 


^ = |i( 1 +cos (f))' if M< 2 / 

I 0, otherwise. 


(2.7) 


( 2 . 8 ) 


The scaling factor A has units of length for forces generated by the ID wall interfaces (for 
which (2.6) approximates a line integral) whereas A has units of area for the 2D biofilm 
region. To ensure that Eq. (2.6) is a consistent representation of the corresponding integrals 
under grid refinement, A is inversely proportional to the number of IB points. Details on the 
precise expression used for A in each case are provided in Section 2.5. 


Step 3: Integrate the Navier-Stokes equations using Chorin's split-step projection scheme: 


a. Compute an intermediate velocity u|y by updating the velocity only for the contribu¬ 
tion from the IB elastic force: 


Af 


= f", 
V 


(2.9) 


(2) (3) 

b. Compute intermediate velocities u' and u). by applying convection and diffusion 
terms using an alternating direction implicit (ADI) approach: 


u (2) -u (1) 

u y u v , <1 » n o n (2) 


K-DXf' I =?D+D- 


Af 


V 


(3) _ (2) 

U ±_ a iL + u n D 0 U (3) | =M D+ D - U (3) 


Af 


V 


( 2 . 10 ) 

( 2 . 11 ) 


The operators D+ and D x refer to the standard first-order forward and backward differ¬ 
ence approximations of the x-derivative, and D x is the standard second-order centered 
difference approximation. Analogous definitions apply for the y-derivative approxima¬ 
tions D+, D y and D°. Equations (2.10) and (2.11) thus represent periodic tridiagonal 
linear systems for the intermediate velocities. 

(3) 

c. Project the intermediate velocity u-. onto the space of divergence-free vector fields by 
first solving a Poisson equation for the pressure p’-j 1 1 

( 2 . 12 ) 

where V;, = (D x ,Dy) is a centered approximation of the gradient operator and the dis¬ 
crete Laplacian • V/ 7 yields a wide finite difference stencil that spans four grid points 
in each direction as described in [50]. Because of the periodic boundary conditions. 
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the pressure Poisson equation is solved most efficiently using a fast Fourier transform 
(FFT). Finally, the velocity projection is completed via the correction 


<j ¥ '='tf ) -jVhP,i- (2-13) 

Step 4: Evolve the immersed boundary to time t n+ \ using 

X' !+1 = X” + Af£>'; +1 4(xy ^ -X$)h x h y . (2.14) 

y 

The IB algorithm explained above is first-order accurate in both time and space. Despite the use of 
second-order differences for spatial derivatives, the spatial accuracy reduces to first order owing 
to the particular choice of interpolation used for spreading the fluid velocity onto the immersed 
boundary [37]. 

This algorithm has the advantage that it is simple and easy to code, although it suffers from 
a fairly strict stability restriction on the time step owing to the explicit treatment of the IB forcing 
term in Step 1 of the algorithm. We implement the above algorithm using MAT LAB®, and when 
that is combined with the extra cost of implementing a realistic detachment criterion (see details 
Section 3.3) our approach is restricted to fairly short-time simulations. The scope of this paper is 
therefore limited to the study of biofilm deformation and detachment in the early stages up until 
a quasi-steady biofilm configuration in reached. Any full-scale implementation of the detachment 
criteria described in this study would therefore benefit from a more efficient IB implementation 
such as the fully implicit approach in [38] or one of the parallel approaches developed in either [24] 
or [62], 

2.4 Discrete representation of walls and biofilm 

In this section, we describe the discretization of the immersed boundaries (walls and biofilm) in 
terms of a network of IB points connected by springs. Section 2.5 will then provide details of the 
force density calculations based on this discrete IB configuration. 

We begin with the horizontal walls at locations y = 0 and H (labeled e-f and g-h in Fig. 1) that 
are each replaced by a ID periodic array of IB points stretching across the domain. Adjacent wall 
points are not explicitly connected to each other, but instead each IB point is connected by a very 
stiff spring to a corresponding tether or target point that initially occupies the same position as the 
IB point. The tether forces generated by these springs constrain the IB points to remain close to 
their target locations and hence mimic a solid wall. The initial IB point spacing h wa u is chosen 
so that h wa ii < \max{h x ,h y ), which helps to control numerical errors that would otherwise lead to 
significant leakage of fluid between IB points [43]. 

Next we consider the discrete representation of the biofilm colony, which is obtained by trian¬ 
gulating the biofilm region, placing IB points at the nodal locations, and constructing a network of 
Hookean springs corresponding to the edges in the triangulation. The spring stiffness is chosen to 
approximate the elastic properties of an actual biofilm material. The biofilm is attached to the lower 
wall by connecting points along the bottom of the colony to a corresponding wall IB point with a 
stiff spring, where each pair of points initially occupies the same location. To ensure that such a 
spring network mimics a mechanically isotropic biofilm, we use an initial triangulation that has 
approximately uniform shape and size of the elements. For this purpose, we use the open-source 
MATLAB package DistMesh of Persson and Strang [41], The average spacing between biofilm IB 
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points is initially chosen to be roughly equal to one-third of the fluid grid spacing. This choice is 
motivated by Vo et al. [57] and ensures that the convective flow inside the biofilm is negligible, 
so that the biofilm acts as an impermeable material. Nevertheless, we note that during the course 
of our biofilm simulations as the biofilm deforms and IB points reach their maximum separation, 
it is possible for portions of the biofilm colony to experience a small apparent permeability. We 
will return to this point when devising a method for approximating the interfacial shear stress in 
Section 3.1. 

Fig. 2 depicts triangulations of two initial biofilm shapes used in this study, one a super-ellipse 
with height Hj, = 75 /mi, and the other a mushroom-shaped biofilm colony with height Hj, ~280pm. 
A coarser mesh is depicted than what is actually used so that the shape and distribution of trian¬ 
gles are evident. IB simulations with a super-ellipse show that the sharp 90 degree corners at 

(a) Super-ellipse shape, with filleted corners at the wall. (b) Mushroom shape. 



280 Jim 



Figure 2: Quasi-uniform triangulations generated by DistMesh for: (a) a super-elliptical biofilm colony with height Hj,=75/(m 
and mean edge length 1.5 f(m; and (b) a taller mushroom-shaped colony with mean edge length 4 Jim. 

the intersection between the biofilm colony periphery and the bottom wall cause localized dis¬ 
continuities in the slope of the biofilm colony-fluid interface to develop at these locations as the 
biofilm colony begins to deform. We believe this is due to the corner-type singularities appearing 
in the fluid that are communicated to the immersed structure via delta function interpolation. In 
actual biofilms, such sharp corners seldom occur and so we instead smooth out the corners using 
a quarter-circular curve or fillet. As depicted in Fig. 2(a), the size of the fillet is kept small so as to 
minimize its influence on either the biofilm or the corresponding elastic forces. 


2.5 Discrete force density calculation 

As indicated earlier, the force density function consists of contributions from two classes of im¬ 
mersed boundaries: horizontal rigid walls (F wa u) and elastic deformable biofilm regions (F/„ 0 ). The 
force contribution in both cases is calculated using a discrete specification that is defined in terms 
of the current configuration of IB points in either walls or biofilm. 
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2.5.1 Top and bottom wall forces 

The discrete wall force arising at any given wall IB point is determined from the stretched state of 
the spring of zero resting length that connects it to the corresponding tether point, yielding a force 
density 


KattJ = *wM (KalU-Kthj)' ( 2 - 15 ) 

where K wa u (g/cm 2 s 2 ) is the spring stiffness, and f and X” f are the coordinates of the wall 
at tether points at time level n. For the stationary bottom wall all tether points are fixed in time, 
whereas the top wall tether points move with a given constant velocity (U wa u, 0) according to 

Kih/ = X ?eth/+(U wa ii,0)At mod (W,H), (2.16) 

where the "modulo" operator ensures that IB points remain inside the domain by imposing the 
periodic boundary condition in the horizontal direction. By choosing K wa u sufficiently large, we 
ensure that the wall IB points do not deviate significantly from their tether point locations through¬ 
out a simulation. For both sets of wall IB points, we set the integral scaling factor in Eq. (2.6) equal 
to the tether point spacing, A = h wa n. 


2.5.2 Biofilm forces 

The force density at an IB point within the biofilm region is determined by summing up the elastic 
force contributions coming from all springs connected to that node in the triangulation. Any given 
spring link is identified by an index pair £,m corresponding to two IB points with coordinates X” 
and X” . The elastic force density contribution at node £ due to spring £,m acts in the direction of 
the vector d" m = X” — X”, joining nodes £ and m. By denoting the resting length of this spring d ( j m , 
the force density F^. , at node £ owing to all springs attached to that node can be written as [1] 


■p/; _ K bio 

r bio,e— ^0 


N b 


m =1 


‘‘L < 


(2.17) 


where jq,, 0 is the spring stiffness coefficient (g/cms 2 ), d° is the average spring resting length and 
d(m = \\df m \\. The symbol ll/ m is a square connectivity matrix of dimension N), x N/, whose entries are 
either 1 or 0 depending on whether or not nodes £ and m are connected. Implicit in this notation 
is the fact that the summation is only done over pairs of nodes that are connected by an edge in 
the network. The scaling factor A in Eq. (2.6) for the biofilm force spreading term is taken equal 
to the average area of a triangle in the biofilm at its rest state, which is equal to the total initial 
area divided by N), (and hence has units of cm 2 ). Finally, following the arguments of Alpkvist 
and Klapper [1], we can ensure that the biofilm deformation is independent of grid refinement by 
scaling the spring stiffness value with the nodal mean distance d° and setting k/, ;o = K^/d 0 . 


3 Interfacial shear stress, drag and lift forces, and detachment 

3.1 Computing forces on the biofilm-fluid interface 

Determining the forces acting on a biofilm colony as it deforms in response to fluid shear is of 
fundamental importance in this study. The most common global measures of hydrodynamic force 
in such FSI simulations are the drag and lift force. In addition, we are interested in calculating the 
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local interfacial shear stress along the biofilm-fluid interface, owing to its essential role in deter¬ 
mining both the mode and extent of biofilm colony detachment. Drag and lift forces are typically 
calculated in the IB framework by summing the corresponding components of the IB force along 
the interface [22,33]. Instead, we follow the approach of Williams, Fauci and Gaver [63] (which we 
refer from this point on as WFG) wherein the traction force is first calculated from the interfacial 
shear stress and then integrated along the biofilm-fluid interface. Our aim in this section is there¬ 
fore to first obtain an expression for interfacial shear stress, and then to derive expressions for the 
drag and lift forces. 

Evaluating interfacial stress in the IB framework is complicated not only by the diffuse nature 
of immersed boundaries owing to regularized IB forces, but also because of the spatial averaging 
inherent in the velocity no-slip boundary condition. If not handled appropriately, both of these 
steps can introduce large errors in interfacial shear stress. This issue was studied by WFG [63] who 
proposed two methods for calculating the tangential interfacial shear stress based on the equation 

tj4n=-pj- (3-1) 

Jump in FS '— 

WS 

Here, X(s,f) is a parametric representation of the biofilm-fluid interface, n = f/If I ^ the unit 
normal vector (directed outward from the biofilm into the surrounding fluid), t is the counter¬ 
clockwise tangent vector, a = — pl + /r(Vu-|-Vu r ) is the fluid stress tensor, and square brackets 
[•] denote the jump in a quantity across the interface. WFG proposed evaluating the tangential 
stress using either side of Eq. (3.1): the left hand side requires calculating the jump in fluid stress 
across the biofilm-fluid interface, and so is referred to as the FS method; whereas the right hand 
side involves local IB force densities, and is called the wall stress or WS method. WFG performed 
IB simulations of 2D Poiseuille flow in a channel, with and without obstructions, and drew the 
following conclusions about the relative merits of these two methods: 

• On a curved boundary, the WS method over-estimates shear stress in comparison with the 
FS method. 

• In order to maximize accuracy with the FS method, the fluid shear stress associated with an 
IB point on the interface should be evaluated at a point located inside the domain, directed 
along the normal vector and separated from the boundary by a distance equal to the fluid 
grid spacing. The fluid grid cell within which this point falls is called the interpolation box 
because it defines a set of 4 fluid grid points that will be used for interpolating the stress. 

• For a curved immersed boundary, the estimate for shear stress at certain IB points can be 
adversely affected when the boundary intersects the interpolation box and a portion of the 
interpolation box lies outside the fluid region. Therefore, only those stresses determined 
using interpolation boxes that lie entirely on one side of the interface should be included, 
and for this purpose WFG designed an exclusion filter that omits any such unwanted stress 
contributions. 

We implemented the WFG exclusion filter approach and found that while it works well for the 
simple immersed boundaries considered in [63], the accuracy of the shear stress calculation de¬ 
grades for the highly curved interfaces that are so common in biofilm applications. We therefore 
developed a modified version of the exclusion filter mentioned above that handles highly curved 
boundaries in a robust fashion. 
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We begin by explaining the WFG exclusion filter in reference to Fig. 3, which depicts the IB 
points X'/' 1 for £=l,2,...,N int lying along a biofilm-fluid interface, ordered in the counter-clockwise 
direction. For each IB point, the outward unit normal m is computed using a local cubic Lagrange 
interpolation procedure as explained in [68]. We then identify two points that will be used to 
evaluate the stress, XE° e ut =X l " t +h x n e and XE" 1 =X l p nt —h x np, which lie on either side of the interface 
located a distance h x along the outward/inward normals from X'" 1 , respectively (assuming here 
that h x = hy). The interpolation boxes or fluid grid cells in which these two evaluation points XE/ U/ 
and XE"' lie are denoted abed and a'b'c'd ', and their respective centroids by XC°/ lt and XC 1 / 1 . Then, 




Figure 3: (a) The exclusion filter uses fluid grid cells abed and a'b'c'd 1 as the two interpolation boxes, determined by 
extending an IB point along the outward and inward normals. For the sake of clarity, IB points lying inside the biofilm are 
not shown here. Distances from the interpolation box to the IB can be calculated in two ways as shown in the inserts: (b) 
measured from the centroid of the interpolation box, or (c) measured as the shortest distance from the edges. 

we calculate the minimum distance from the biofilm-fluid interface to the centroids of the two 
interpolation boxes (referring to Fig. 3(b)): 

MD° ut = min () and MDf = min (d^), (3.2) 

where m ranges over 1,2and Df" and D' f n m represent the distance between IB point m and 
the centroid of the two interpolation boxes associated with the 1 th IB point. 

The main principle behind the WFG exclusion filter [63] is to base stress calculations only on 
those points that are most representative of the fluid flow around the biofilm colony, first by ex¬ 
cluding interpolation boxes that are cut by the interface, and second by choosing stress evaluation 
points XE that are as far from the interface as possible (so that they are least affected by the regu¬ 
larized IB force). To this end, we identify IB points i for which MD P has a local maximum, while 
also requiring that no interpolation box be used for more than one stress calculation. 

As mentioned above, when this WFG filter is applied to biofilms with highly curved interfaces, 
the accuracy of the stress approximation degrades. One cause of this degradation is that certain 
portions of the interface may end up with few IB points included, leading to poor resolution. To 
address this problem, we propose a modification of the WFG exclusion filter that is implemented in 
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Algorithm 1 for stresses on the outer (fluid) side of the interface; only minor changes are required 
for stresses on the inner (biofilm) side. This new filter makes the following three modifications to 
the basic WFG exclusion filter: 

• In step 8, check whether each interpolation box is cut by the biofilm-fluid interface. 

• In step 11, compute the minimum distance between the biofilm-fluid interface and the edges 
of the interpolation box abed as shown in Fig. 3(c), rather than the centroid in the original 
filter. This provides additional resolution by differentiating between nearby IB points that are 
equidistant from the centroid, since points are not equally-spaced from the edges of the box 
(compare Figs. 3(b,c)). The minimum distance calculation from the edge of an interpolation 
box is performed using a fast and elegant algorithm from [4] that introduces no significant 
additional cost. 

• In step 18, the strong local maximum criterion is relaxed and we instead introduce a new 
user-specified filtering parameter e min . This change ensures (for an appropriate choice of 
e m i„) that we retain some IB points that would otherwise have been rejected by the original 
WFG filter. 


Algorithm 1 Modified WFG exclusion filter, which determines a list of IB points (pt_list) and 
interpolation boxes (box_list) for stress calculations. 

1: for all interface IB points X'" ! do 
2: Compute the outward unit normal rq. 

3: Locate the evaluation point XE"" f = X(' if + /i r n/. 

4: Identify the interpolation box abed corresponding to the fluid cell containing XE(’“ / . 

5: Add abed to the sorted list box_list and remove duplicates. 

6: end for 

7: for all interpolation boxes abed G box_list do 

8: if abed is intersected by any of the N mt — 1 line segments comprising the biofilm-fluid inter¬ 

face then 

9: Remove abed from box_list. 

10: else 

11: Determine the minimum distance MD° ut from the edges of box abed to the biofilm-fluid 

interface, and the corresponding IB point i. 

12: Store MD° ut and abed in the data structure for IB point £. 

13: end if 

14: end for 

15: Initialize pt_list to contain a list of all interface IB points. 

16: for all Xg G pt_list do 

17: Let abed be the interpolation box corresponding to Xg. 

18: if either of the neighbouring IB points £±1 has the same interpolation box abed 

or the minimum distance MDg< e, mn then 
19: Remove X^ from pt_li st. 

20: end if 

21: end for 


After executing Algorithm 1, we obtain a list of IB points at which fluid stress tensor compo¬ 
nents can be estimated inside corresponding interpolation boxes. We first determine approxima- 
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tions of the velocity derivatives at the corners of each interpolation box using one-sided second- 
order differences, where the points included in the difference stencils are chosen to lie entirely 
inside (or outside) the biofilm colony as determined by the unit normal to the interface. We then 
take the velocity derivatives and pressures, and apply bicubic interpolation [45] to determine cor¬ 
responding values at the stress evaluation point, which are then combined to obtain the stress 
tensor components, cru, and C \2 = (J 2 \ . In most cases, the set of IB points selected for computing 
the stress tensor components inside and outside the biofilm colony are different; therefore, com¬ 
puting stress jumps requires that stresses be interpolated onto points in each set. We may then 
obtain all necessary values of the the stress jump [cr] and the tangential component of the interfa¬ 
cial shear stress, t- [cr] n. Finally, the drag and lift forces are computed by integrating the jump in 
shear stress along the biofilm-fluid interface using 

©=/ ( ^ 2) -(|Sl m) ds - (3 - 3) 

3.2 Computing the averaged equivalent continuum stress inside the biofilm 

In our IB model, the biofilm continuum is replaced by a network of discrete springs wherein the 
elastic restoring forces arising from stretched /compressed springs take the place of stress and 
strain in an real elastic continuum. The most natural way to simulate biofilm detachment within 
such a spring network representation is to cut any spring links for which the local strain exceeds 
a critical value, as explained in [1,27]. However, this approach suffers from several drawbacks. 
First of all, the force resulting from stretching or contraction of a ID spring element cannot accu¬ 
rately capture the actual strain in an elastic continuum and consequently there is no direct way 
to determine a critical spring strain threshold based on measured biofilm mechanical properties. 
This contrasts markedly with other approaches such as [7,13,44] that discretize the solid mechan¬ 
ics equations directly and employ a more physically realistic von Mises yield stress criterion to 
initiate detachment. In addition, there is no reliable way in the spring network approach to de¬ 
termine spring parameters so as to ensure that different triangulations exhibit similar detachment 
dynamics under the same flow conditions. 

In order to bridge this gap between continuum mechanics-based biofilm models and the dis¬ 
crete IB spring-based model, we introduce the notion of a stress tensor defined at each node in the 
network. This is accomplished by assuming that there is an equivalent continuum representative 
elementary area or REA surrounding each node, within which we compute an average value of the 
stress tensor components in terms of the spring forces acting on that node. The primary motivation 
for this definition comes from the Discrete Element Method (DEM) for computing microstructural 
stress in a granular medium [5,20]. In the DEM, the dynamics of a granular medium are deter¬ 
mined by treating each grain separately and solving the governing force balance equations under 
the combined action of grain contact forces, body forces and external forces. In place of grains we 
have IB points, and our Hookean spring connections replace the contact forces between the grains. 

3.2.1 Constructing the REA around each IB point 

In contrast with DEM simulations of granular media that identify an REA with a Voronoi cell 
constructed from a Delaunay triangulation [5], we employ instead a control volume finite el¬ 
ement method construction [59] that is based on a (non-Delaunay) triangulation generated by 
DistMesh [41]. For each biofilm point labelled I in the triangulation, the corresponding REA is 
constructed by joining with straight lines the centroids of all surrounding triangles to the mid¬ 
points of the corresponding edges (springs) emanating from node I as depicted in Fig. 4. If there 
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are m springs connected to node I, then the corresponding REA is a polygon with 2m sides. By this 
construction, we note that the area A of the REA around a node (which is required for calculat¬ 
ing the stress tensor later in Eq. (3.12)) is one-third of the total area of all triangles surrounding the 
node. This REA construction extends in a straightforward manner to volumes in three dimensions. 



Figure 4: The representative elementary area (REA) surrounding IB node I corresponds to the polygonal region V (shown 
using dashed lines) connecting the centroids ( a,b,c,d,e ) of the surrounding triangles to mid-points {j,k,l,m,n) of the edges 
emanating from node I. 


3.2.2 Computing the equivalent continuum stress tensor 

Consider the REA surrounding a point X 1 in the triangulation, pictured as a dashed polygonal 
region in Fig. 4. Denote the REA and its boundary by V and d'P respectively, and let A refer to 
the area of V. Suppose that an equivalent continuum material occupies the REA, with stress field 
having Cartesian components for i,j= 1,2. Then the average stress over V is 




which can be manipulated to obtain 


[K? x ‘^~ x ‘ <%* 


dA, 


(3.4) 


where the subscript ", k" denotes a /c-component derivative and the Einstein summation conven¬ 
tion is assumed for repeated indices. The divergence theorem may then be applied to the first term 
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to get 

^ (£, 1 °kj X i U k dS ~ fj v X i a kj,k dA ) • ( 3 ‘ 5 ) 

Now consider the two types of force that can act on the REA: a surface traction force T)(x) that 
acts at points on the REA boundary, and a body force gj(x) acting at interior points. Imposing a 
force balance on the boundary yields 

o^m(x) = Tj(x) (3.6) 

at points x G dV, where n,- denotes the outward-pointing unit normal vector to dV. We note here 
that the convention in solid mechanics is to use the inward normal (which ensures compressive 
stresses are positive), however we break this convention for the sake of consistency with rest of the 
text. Balancing forces in the interior of the REA gives 

< 7 ij,i+ Pgj=P a j' ( 3 - 7 ) 

where Uj is the /-component of acceleration. Substituting Eqs. (3.6)-(3.7) into (3.5) then yields 

= \{ir P XlT ’ dS + JJ v P Xi ^i- a i) dA ) ■ (3-8) 

We now introduce notation for IB nodes X“, a = J,K,L,M,N, that are immediate neighbors 
of node I in the triangulation shown in Fig. 4, along with corresponding edge vectors E^ for /3 = 
j,k,l,m,n directed outward along springs (with E^=X^— X 1 when jS=/, for example). The edge mid¬ 
points are then denoted by Z-’ = X I + E^/2, with corresponding spring forces F ,,; . If we associate 
the boundary traction force T(x) for x G dV with the spring forces F^, then the boundary integral 
term in (3.8) may be rewritten in component form as 



TjdS = 


E z h(- 

pedV 


(3.9) 


Substituting this expression into (3.8) yields 


cf.? = - V 
v A ' 


K pedV 


yH 7Tr 

z < t i 


- fj[pP x i{gj-a j) dA 


(3.10) 


We assume in this paper that the biofilm is neutrally buoyant and that the equivalent con¬ 
tinuum stress is computed in a steady-state configuration for which inertial forces are negligible; 
consequently, the integral term in Eq. (3.10) is zero. The remaining summation term can be further 
simplified by substituting = X z + E^/2 and using the equilibrium condition 

Xf E F f = O' (3.11) 

pedV 


to obtain 



2 A 


pedp 


F r 


(3.12) 
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We note that this averaged equivalent continuum stress tensor is guaranteed to be symmetric (u^= 

17j .?), which should be contrasted with the analogous derivation for granular media where the DEM 
approach leads to a non-symmetric stress tensor owing to contact forces with a nonzero moment 
about grain centers [5]. This does not happen in our IB spring network because elastic spring forces 
always act along lines connecting IB points and hence do not generate any such moments. 

For our choice of polygonal REA around node I, the area A in (3.12) is equal to one-third of 
the total area of all triangles surrounding the node. If we apply our stress calculation method to 
a biofilm that only deforms and experiences no detachment, the number of springs connected to 
each node is constant and a simple data structure can be used to store the information needed 
to compute stress. Furthermore, the area A and edge vector components E? are easily updated 
using the current IB node coordinates. However, in the more complicated case with detachment, 
extra geometric information must be stored along with the IB point data; for example, we must 
keep track of the changing number springs connecting each IB node as well as the number of 
active triangles around each node. The necessary changes to the code and data structures are 
straightforward, and the added computational cost is negligible in comparison to that for the IB 
algorithm. 

3.3 Biofilm detachment criterion 

As mentioned earlier, we handle biofilm detachment using a yield stress criterion from the von Mises 
stress theory, which is one of many approaches used to model failure of a ductile material subjected 
to external loading [3]. The von Mises yield stress criterion has been employed successfully in the 
context of biofilm modeling [7,13,44], even though the biofilm composite (composed of cells, EPS 
and fluid) is not strictly a ductile material. In two dimensions, the von Mises yield stress at any 
IB point inside the biofilm region can be expressed in terms of the averaged equivalent continuum 
stress tensor components in Eq. (3.12) as 

(<7von) 2 = (o : l!) 2 -^22 + (^22) 2 + 3 (^12) 2 • ( 3 • 13 ) 

The value of cr von is then compared to some threshold yield stress, which is a measure of the biofilm 
cohesive strength. If a von exceeds the threshold, then detachment is initiated by severing all springs 
connecting it to neighbouring IB points. This should be contrasted with other approaches based 
on edge strain [1,27], wherein springs are cut individually based on some threshold strain value 
and an IB point is only detached from the colony when all springs connected to it are severed. 

The detachment process in the case of biofilms is complicated somewhat by the fact that yield 
strength varies throughout the biofilm colony. For example, the strength with which the base of 
the biofilm colony adheres to the substratum, which we call the biofilm adhesive strength E,,,//,, is 
several orders of magnitude larger than the cohesive strength of the bulk biofilm material. Further¬ 
more, the bulk cohesive strength also varies since the portion of the colony nearest the biofilm-fluid 
interface has a stress threshold that is significantly smaller than the value Z'"^ in the interior. 
Consequently, we have three threshold values satisfying E^ < E™^ < Z which leads to a natural 
separation of the biofilm into three zones - substratum, interface and interior. 

With this in mind, we propose the following detachment strategy. First, for any given IB point 
X/ we calculate the shortest distance from point £ to the substratum and to the biofilm-fluid inter¬ 
face, denoted by Df lb and O l ’ f xt respectively. We then select a yield stress criterion to be imposed 
by determining which zone the IB point belongs to: 
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Zone 1: consists of all IB points near the substratum that satisfy Df lb ^ e sub , where e sub is a user- 
specified parameter. In this case, the point detaches whenever a von 'L ndh . 

Zone 2: consists of any of remaining IB points near the biofilm-fluid interface that satisfy D‘’ xt ^ 
e ext . Here, the point detaches when a von ^ 

Zone 3: consists of all remaining interior biofilm points, which detach if cr von ^ . 

Calculating the distance Df lb is trivial because the bottom wall in our numerical simulations is 
parallel to the x-axis. However, the calculation of D exi is more involved owing to the irregular 
biofilm shape and also because the colony deforms in time, hence requiring that D‘j xl be recalcu¬ 
lated in each time step. Therefore, care must be taken in order to design an efficient algorithm for 
estimating D ext ; for this purpose we employ a finite element-based signed distance function for 
triangles developed in [17], which is an extension of the fast marching method. The cost of this 
algorithm can be optimized by only calculating the distance function at IB points lying within a 
narrow band near the interface. As an illustration, Fig. 5 shows sample contour plots of D‘j xl for 
two different biofilm colony shapes. 



(a) (b) 

Figure 5: Contours of the distance function D| If (labels in cm) from any given IB point to the biofilm-fluid interface, 
computed using the algorithm in [17], Two biofilm colonies are shown near steady-state: (a) the elliptical configuration 
SUP75, and (b) a mushroom-shaped colony. 


The treatment of biofilm detachment is performed after completing Step 1 (the force calculation 
step) in the immersed boundary algorithm. The steps in the yield stress based detachment process 
are summarized in Algorithm 2. We note that in the implementation outlined here, we make use 
of three integer arrays of "status flags" - named *_STATUS with * = INODE, ISPRING, ITRI - one 
each for IB nodes, springs and triangles. These flags are set to either 0 or 1 depending on whether 
the status is detached or active, respectively. 















21 


Algorithm 2 Biofilm detachment algorithm based on equivalent continuum stress and the 
von Mises yield stress criterion. 

1: For each IB node £, compute the distance from the substratum, D s " b . 

2: Identify IB nodes i that lie on the interface, and assign D e „ xt = 0 there. 

3: At all remaining IB nodes, compute the distance function D' ; f xt using the algorithm in [17]. 

4: for all IB nodes X/ that are active (with INODE_STATUS = 1) do 

5: Compute the REA area A as the sum of areas of all active triangles (with ITRI_STATUS = 1) 

neighbouring node £. 

6: Use all active springs (with ISPRING_STATUS = 1) to calculate the equivalent continuum 

stress components erf? from Eq. (3.12). 

7: Compute the von Mises yield stress a von using Eq. (3.13). 

8: if ^ e sub then 

9: Let THRESHOLD = IL adh . 

10: else if D ext ^ £ ext then 

11: Let THRESHOLD = 

12: else 

13: Let THRESHOLD = E^. 

14: end if 

15: if lT von ^ THRESHOLD then 

16: Let IN0DE_STATUS = 0 (detached). 

17: For each spring adjacent to this node, set ISPRING_STATUS = 0. 

18: Use IN0DE_STATUS and ISPRING_STATUS to update the status of the triangle associated with 

the detached node to inactive (ITRI_STATUS = 0). 

19: end if 

20: end for 
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4 Numerical simulations 

4.1 Model parameters 

Table 1 summarizes all parameter values used in our numerical simulations of biofilm-fluid inter¬ 
action. The parameters are separated naturally into the following categories: 


Table 1: Parameter values for the various numerical test cases. 


Description 

Values 

Fluid domain: 

Domain height 

H = 3xH b 

Colony spacing 

D b = 50, 150, 250, 400 pm (W b ,3W b ,5W b ,8W b ) 

Biofilm shape (width and height): 

Semi-circle: 

W b = 40 pm, H b = 20 pm (SEMI20) 

Super-ellipse: 

W b = 50 pm, H b = 25 pm (Sup25) 

W b = 50 pm, H b = 50 pm (Sup50) 

W/, = 50 pm, H b = 75 pm (SUP75) 

Fluid/biofilm grid: 

Fluid grid spacings 

h x = h y = 0.5 pm (Semi20, Sup25) 
h x =hy = 0.75 pm (Sup50, Sup75) 

IB wall point spacing 

h waU = \mm(h x ,h y ) 

Biofilm spring rest-length 

d® = 0.15 pm (Semi20, Sup25) 
d° = 0.225pm (Sup50, Sup75) 

Fluid/biofilm material properties: 

Fluid density 

p = 1.0 g/cm 3 

Fluid viscosity 

^ = 0.01g/cms 

Shear rate 

G = 0.625 s- 1 

Biofilm spring stiffness 

K b io = 0.75d°, 7.5d°, 75d° g/cm 2 s 2 

Wall spring stiffness 

K waii = 10 5 g/cm 2 s 2 


• Biofilm geometry: We took four different biofilm colony shapes, one a semi-circle of radius 
20 pm and width IV/, = 40 (labeled SEMI20), and three (semi-)super-ellipses having width 
Wf, = 50/(m and height H b = 25, 50 and 75}im (labeled Sup25, Sup50, Sup75 respectively). 
These dimensions are representative of typical biofilm colonies and also capture a range of 
aspect ratios observed in early stages of biofilm colony growth. 

• Fluid domain: The vertical spacing H between top and bottom channel walls is set to triple the 
height of the biofilm colony (H = 3 H b ) in order to minimize boundary effects. Simulations 
reveal that at our results are relatively insensitive to changes in H. The width of the fluid 
domain is W b + D b , and values of O b = W h , 31V/,, 5IV/, and 81V/, were used to study the effect 
of colony spacing. 

• Fluid grid: We chose relatively small values of fluid grid spacing h x and h y that permit ac¬ 
curate resolution of the biofilm colony. In particular, we aimed to ensure that recirculating 
eddies arising from flow separation are well captured. A constant time step of A t = 10 3 s 
was used in all simulations. 
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• Biofilm grid: The mean spacing d° between biofilm IB points was chosen to satisfy d° ^ 
^ min (h x , h y ) as in [57] so as to avoid numerical errors due to leakage of fluid between IB 
points. This value of d° is provided to the DistMesh code as a measure of average edge length 
for the biofilm triangulation; for example, the SUP75 biofilm colony with d° = 0.225 pm yields 
a triangulation with 81,449 IB nodes and 222,659 edges. 

• Fluid material properties: We take fluid parameters consistent with water: p — 1 g/cm 3 and 
p = 0.01 g/cms. The shear rate was set to G = 0.625 s 1 for all simulations, which is high 
enough to induce large deformations in weak biofilm colonies while also attaining a steady 
state over a relatively short time period (roughly 2-8 s). This shear rate corresponds to a 
Reynolds number Re=pGWfi/p = 1.5625x 10 -3 for the case W/, = 50ym, where we have used 
the width of the biofilm colony as a length scale. 

• Biofilm material properties: To mimic weak biofilms with varying mechanical strength, we 
choose several values of the IB spring stiffness corresponding to k^- o =0.75, 7.5 and 75, where 
we recall that k&; 0 = jq,; 0 / d°. The wall spring stiffness K wa u = 10 5 is chosen much larger so that 
wall points do not move appreciably. 

In summary, we consider four different biofilm shapes, four colony spacings and three values of 
the spring constant, corresponding to a total of 48 simulations with a single value of shear rate 
G = 0.625 s- 1 . 

4.2 Model validation: Channel with a rigid bump 

Before applying our immersed boundary algorithm to the biofilm test cases outlined in the previ¬ 
ous section, we first validate our numerical approach using a simpler set-up consisting of a rectan¬ 
gular channel containing a rigid semi-circular bump on the bottom wall. The same problem was 
considered in Williams et al. [63] as an illustration of their exclusion filter approach for interfacial 
stress calculation. 

The problem geometry is shown in Fig. 6, consisting of a channel of width H constructed from 
two parallel horizontal walls along y = y and y = The channel is embedded within a larger 
rectangular fluid domain of size L x 2 H and periodic boundary conditions are imposed the outer 
boundaries. A constant fluid body force fg(x) = (A P/L, 0) is applied at each fluid grid point 
inside the channel, which in the absence of any obstruction would generate a parabolic Poiseuille 
flow with pressure difference A P across the channel from left to right. However, an obstruction is 
introduced within the channel consisting of an immersed boundary in the shape of a semi-circle 
of radius yy and centered at (y,y + |). The corners connecting the bump to the bottom wall are 
smoothed using quarter-circular fillets with radius as shown in Fig. 6(b) that serve to regularize 
the IB shape and avoid flow irregularities near sharp corners. 

The channel walls and semi-circular obstruction are represented using IB points, each of which 
is connected by a single spring to a tether point that is fixed in space, and the spring stiffness 
is chosen large enough that the IB points do not move appreciably. The spacing h wa u between 
adjacent IB points is chosen so that h wa u ^ |min (h x , hy), which aims to minimize any numerical 
errors arising from leakage of fluid across the immersed boundaries. 

We choose parameter values the same as in [63], namely H = 0.1 m, L = 0.8 m, R = 0.05 m, 
A P/L = — 10 5 kg/m 2 , // = 25kg/ms, p = lkg/m 3 and k = 10 9 kg/m 2 s 2 . In contrast with the CGS 
units used in the rest of this paper, we employ MKS units in this section only for ease of comparison 
with the results in [63]. The fluid domain is discretized on a uniform grid of 512 x 128 points and 
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Figure 6: (a) Computational domain for the channel flow with a rigid bump. The fluid domain is a box of size Lx2H 
containing two parallel, horizontal walls at height y= ^ and ^r-. Attached to the bottom wall is a filleted (smoothed) 
semi-circular bump. The flow is driven from left to right through the central channel by imposing a constant (positive) 
fluid body force, (b) Zoomed view of the semi-circular obstruction. 

















25 


we use a constant time step A t = 0.16 ph^ n / pi. The Reynolds number based on channel height is 
Re = pAPH 3 / (8/r L) ~ 0.02, which indicates that inertial effects are negligible and permits us to 
make comparisons with the Stokes flow solution from Gaver and Kute [21]. 




(a) 


(b) 


Figure 7: (a) Dimensionless interfacial shear stress computed along a semi-circular obstruction in a channel. For comparison, 
results are included from another IB approach [63] and the boundary element method [21]. (b) Close-up view near the apex 
of the bump where shear stress attains a maximum. 


The computed interfacial shear stress, non-dimensionalized by A P, is shown in Fig. 7(a) along 
the central portion of the bottom wall including the semi-circular bump. Results are presented for 
two values of the exclusion filter parameter, e m i n =0.15h x and 0.4/z x . Also included in this figure are 
numerical results computed with two other methods, namely the WFG method with an FS-based 
exclusion filter [63], and the boundary element computations of Gaver and Kute [21]. The results 
from our revised exclusion filter are within approximately 5% of WFG's results, with the differ¬ 
ences most pronounced on either side of the crest of the bump (see Fig. 7(b)). The correspondence 
with WFG's results improves as increases from 0.1 5h x to 0.4 h x , which can be explained as 
follows: even though the number of IB points retained by the filter decreases with increasing e m {„, 
the quality of those points is high because they are separated further from the smearing effects of 
the interface. 

A non-dimensional flow rate can be computed by integrating the computed velocity vertically 
along the channel inlet 
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12/zL 
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(4.1) 


yielding a value of Q = 0.554 that agrees exactly with the result of WFG [63] for the same grid 
resolution. We also compute the dimensionless drag 
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where Ti is the x-component of the traction force and F represents the portion of the bottom wall 
corresponding to the semi-circular bump. Our simulations yield /£, =4.662 and 5.195 for filter 
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parameters e mm =0.15h x and 0.40/z x respectively, which should be compared with the WFG result of 
5.617. This discrepancy of 7-17% is acceptable in view of the increased flexibility we gain from our 
modified filter in terms of being able to compute stress along strongly-curved biofilm interfaces. 


4.3 Simulating biofilm deformation: Flow structure and forces 

In this section, we investigate the response of deformable biofilm colonies to a shear flow by study¬ 
ing the effect of changes in various parameters on the biofilm shape, flow structure, hydrodynamic 
drag and interfacial shear stress. 


Undeformed 
D b = 50 
D b = 150 (xm 
D b = 250 
D b = 400 
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D b = 50 |xm 

■ D b = 150 (xm 
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Figure 8: Initial (undeformed) and steady-state (deformed) biofilm shapes with = 0.75 as the spacing between colonies 
Dj, is varied, (a) Test cases SEMI20 and SUP25. (b) Test cases SUP50 and SUP75. 

We begin by varying the inter-colony spacing 1% for the four initial colony shapes SEMI20, 
SUP25, SUP50, SUP75. Fig. 8 depicts the initial and final (steady-state) biofilm profiles for values 
of D h lying between 50 and 400 /xm, holding the stiffness parameter K^io = 0.75. The extent of 
the deformation clearly increases as the biofilm height is increased, with the largest deformations 
occurring for the SUP75 colony. This behaviour is physically reasonable because longer structures 
are not only more flexible but also extend further into the shear flow where they experience a 
higher flow velocity. The extent of deformation also increases with the spacing parameter D/„ 
which is to be expected since the shear flow is more able to impinge between colonies having a 
greater separation. 

Note that our simulated biofilm colonies appear to simply shear to the right without exhibit¬ 
ing any of the elongation or vertical lifting that is observed in the numerical simulations of Vo 
et al. [57,58]. This discrepancy can be attributed to several sources: (a) Vo et al. used different 
colony shapes having a sinusoidal profile that is much taller (corresponding to H = 180 f/m); (b) 
their spring stiffness is roughly 100 times larger than ours; and (c) they also considered a much 
faster flow corresponding to Re = 230 (based on hydraulic diameter of the square capillary reactor 
as length scale) whereas we have Re = pGW^ /p ~ 10 3 (based on a much shorter length scale Wj, 
corresponding to the biofilm width). It was also suggested in [58] that a significant contributor 
to elongation from lifting at high enough Reynolds number is inertial effects arising from biofilm 
deformation as well as the nature of the flow surrounding the colony. 
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Figure 9: Final (steady state) biofilm shapes for different values of Kj ! 0 and constant colony spacing Dj, = 250/;m. The 
super-ellipse test cases SUP25, SUP50 and SUP75 are shown. 


We next investigate the effect of changes in the spring stiffness on the steady-state deformation, 
taking values of Kbio = 0.75, 7.5 and 75 for three different initial colony shapes while holding D/, = 
250 fim constant. Fig. 9 depicts the final deformed shapes from which we observe that at the 
lowest shear rate, even a relatively modest value of biofilm stiffness Kbio ~ 7.5 is sufficient to resist 
deformation. Indeed, it is only when stiffness is reduced to ic^o = 0-75 that any significant bending 
of the biofilm colony occurs. 

The drag force acting on the biofilm at steady state is then computed for all of simulations 
above. Fig. 10 plots the drag force in each case as a function of colony spacing D&. Here, the 
drag force has been non-dimensionalized using fp =/d / (f/GWj,), where the reference value JiGWb 
can be thought of as the force exerted by a linear shear flow with shear rate G acting on a very 
thin biofilm colony with H/, -C Wj,. The corresponding lift force is not shown because lift is much 
smaller than drag (by at least a factor of 10), not to mention that lift force is much less affected 
by changes in colony spacing. We observe from Fig. 10 that /5 is an increasing function of D/„ 
with the rate of increase being largest for the weakest biofilms having Kbio = 0.75. In particular, 
as Db increases from 50 to 400 }im, the drag force increases by a factor of 50% for cases SUP25 
and SUP50, with the largest colony in SUP75 exhibiting a roughly 100% increase in drag. These 
results should be contrasted with the simulations of rigid biofilms in [52] where the drag force was 
non-monotonic, attaining a local minimum at some intermediate value of Db . Because the drag 
fails to level out at the upper end of the range Wj, A D/, ^ 8W/, considered here, we would have to 
explore significantly larger values of colony separation in order to obtain results consistent with 
an isolated biofilm colony. Because a much larger domain size would be required, we have not 
investigated this high-Dj, regime for reasons of high computational cost. 

Another observation is that for stiffer biofilms (Kbio = 7-5 and 75) organized into the closest- 
spaced colonies (Dj, =50/;m), the drag force increases by at most 20% as the colony height increases 
from SEMI20 to SUP75, whereas the drag on the weaker biofilms (Kbio = 0.75) increases by more 
than 50%. This should be contrasted with the widest-spaced colonies (D;,=400/rm) where the drag 
increase for stiff biofilms is roughly 100% and over 150% for weak biofilms. In the latter case, the 
majority of the drag increase occurs when the biofilm height increases from 50 to 75 }im, which 
cannot be accounted for by the increase in colony surface area alone (i.e., perimeter in 2D). We 
therefore conclude that biofilm colonies may be able to grow into tall structures even if they are 
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Figure 10: Variation of dimensionless drag force (at steady state) with colony spacing Dj. Curves are shown for four 
different initial colony shapes and for Kj !O = 0.75, 75. 

weak mechanically because of the protection they gain from having other colonies in close spatial 
proximity 

More insight into the causes and extent of fluid-induced deformation can be derived by visual¬ 
izing the flow structure using path-lines or fluid particle trajectories. Fig. 11 depicts four path-line 
plots for the SUP75 case, corresponding to biofilms that are weak/stiff (with k/, ;o = 0.75, 75) and 
spaced closely/widely (D^ = 50, 400). For the stiffer biofilm, the flow over the widest-spaced 
colony in Fig. 11(d) clearly exhibits flow separation both upstream and downstream of the colony. 
As the spacing between the stiff biofilm colonies is reduced, the up/downstream eddies merge to 
form a single large eddy as pictured in Fig. 11(c), which is similar to what has been observed in 
numerical simulations of rigid biofilms at high shear rates [52,66]. 

In the case of a weak biofilm, the eddy structure and dynamics are more interesting. In narrow¬ 
spaced colonies (D/, =50) the deflection of the weak biofilm causes the single eddy located between 
the neighbouring colonies to become distorted and "climb" the upstream face as seen in Fig. 11(a). 
For more widely-spaced colonies (D b = 400) with two distinct eddies present, the increased defor¬ 
mation in the weak biofilm causes the upstream eddy to shrink while the downstream eddy grows, 
as shown in Fig. 11(b). 

We remark that for the weak biofilms in Figs. ll(a,b), path-lines clearly traverse the interior of 
the colony which corresponds to a very slow flow (slower by several orders of magnitude than the 
flow outside the biofilm region). The spring network making up our simulated biofilm therefore 
behaves like a porous medium, which has a very small permeability [51] that can be attributed to 
small volume conservation errors that are well-studied in the context of the IB method [23]. The 
porous flow is so small that it has minimal effect on the biofilm deformation or the flow, but it 
cannot be ignored when computing shear stress along the interface using the FS method described 
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Figure 11: Pathlines for SUP75 showing the flow structure and biofilm deformation as a function of stiffness and colony 
separation. 
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earlier in Section 3.1. The porosity is particularly important along portions of the biofilm-fluid 
interface that experience flow separation adjacent to recirculating eddies on the upstream and 
downstream faces of the colony. The effect of the weak porous flow inside the biofilm colony is 
accounted for in the FS method by including its contribution to the jump term in Eq. 3.1. This 
requires performing the filtering operation twice: once to select IB points suitable for calculating 
the fluid stress component outside the biofilm colony, and second time for the stress inside. As 
mentioned earlier in Section 3.1, different IB points are selected for the stress calculations inside 
and outside the biofilm colony, necessitating an interpolation between the two sets of IB points. 

This is in contrast with the stress calculations of WFG [63], where the interior flow was assumed to 
be weak and its contribution to the interfacial stress was neglected. This assumption did not affect 
their calculated interfacial shear stress, since there was no flow separation at the low Reynolds 
numbers they considered. 

To conclude this section, we investigate the interfacial shear stress distribution along the biofilm- 
fluid interface, which is depicted in Fig. 12 for weak/stiff biofilms (k/„ o =0.75, 75) that are closely/widely 
spaced (Dj, = 50 and 400 pm). The shear stress is plotted against the IB point index numbered from 
left to right along the interface (note that the total number of IB points increases with the size of 
the colony). In all cases, the shear stress is non-dimensionalized using the reference value pG, 
which corresponds to the steady uniform shear flow that would occur in the absence of any ob¬ 
stacle. We adopt a sign convention that assumes stress is positive when the flow adjacent to the 
biofilm-fluid interface is in the same direction as the primary channel flow (i.e., from left to right). 

A zero shear stress indicates a point where flow separates and a recirculating eddy attaches to the 
biofilm surface. 

For each test case depicted in Figs. 12(a)-(d) there are four curves, which can be separated into 
two pairs having similar shape: one pair corresponding to weak biofilms with high shear stress, 
and the second to stiff biofilms with low shear stress. Within each pair of curves, increasing the 
colony spacing from 50 to 400 pm causes the maximum shear stress to increase but leaves the 
general shape of the stress curve unchanged; however, there is a slight downstream shift of the 
location of the maximum stress for the three super-ellipse cases SupNN. 

In summary, these results indicate that both the magnitude of the shear stress and its variation 
along the interface can change significantly when the biofilm colony deforms. We will see next 
that this has important implications for biofilm detachment. 

4.4 Simulating biofilm detachment using equivalent continuum stress 

To demonstrate the effectiveness of our equivalent continuum stress-based detachment strategy, 
we now consider a different colony shape pictured in Fig. 1 . This mushroom-shaped colony fea¬ 
tures wide head and base sections connected by a relatively narrow stem. The reason for using 
this shape is two-fold: first, these long and thin structures are more realistic especially during the 
advanced stages of deformation just before a detachment event; and second, that it allows com¬ 
parison with other IB studies of biofilm deformation and detachment that use similar mushroom 
shapes [1,27]. Based on solid mechanics principles and experimental evidence, we know that when 
such an elongated colony is subjected to a sufficiently strong shear flow, it will rupture at a loca¬ 
tion inside the narrow stem region where the cross-sectional area is smallest. Consequently, this 
scenario serves as a simple validation of our detachment algorithm. 

The precise mushroom shape used in our simulations is extracted digitally from [1, Fig. 3a] 
which in turn was taken from experimental images. We use the following values of the physical 
and geometric parameters: H = 300 pm, Dj, = 240 pm, h x = h y = 5 pm, d° = 1.25 pm, G = 0.625 s 1 
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k bio /d 0 = 0.75 (D b = 50pm) 



(a) Case Semi20 


k blo /d 0 =0.75 (D b = 50 |tm) 



(b) Case Sup25 



k bio /d 0 = 0.75 (D b = 50 pm) 



(c) Case Sup 50 (d) Case Sup75 

Figure 12: Dimensionless interfacial shear stress plotted as a function of IB point index, numbered from left to right along 
the biofilm-fluid interface. Results are shown for cases SEMI20, SUP25, SUP50 and SUP75, with two values of colony 
spacing (D;, = 50, 400//m) and two stiffness values (?C{,; 0 =0.75, 75). 
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and Kfo'o = 15. As before, DistMesh is used to triangulate the initial biofilm region, yielding a total 
of 6155 IB nodes connected by 6500 edges (springs). In the remainder of this section, we present 
results using two different detachment strategies: the first based on spring strain; and the second 
based on equivalent averaged continuum stress. In particular, we aim to identify the drawbacks 
of the strain-based approach that in turn highlight the advantages of our new equivalent averaged 
continuum stress approach developed in Section 3. 

4.4.1 Spring strain-based detachment 

We start by considering a detachment strategy based on spring strain, = df m / d° £m — 1, where d f m 
represents the length of a spring and d° (m is the corresponding unstressed (or resting) length. The 
central parameter in this detachment model is the critical strain, £ max , beyond which the spring 
will break. Other studies employing a similar detachment criterion [1,27] have used £ max = 1, 
which coincides with detachment occuring when a spring is stretched to twice its resting length. 
The reason for this choice of critical strain was not justified, even though it is evident on physical 
grounds that £ ma x should not be constant but rather depend on the strength of the biofilm matrix 
(which in our IB model is expressed by K},i 0 ). 

As an illustration of strain-based detachment, we simulate flow over the mushroom-shaped 
biofilm using the parameters indicated above. Fig. 13(a) depicts the resulting biofilm configuration 
at four equally-spaced times between t = 0 and 1 s. Unlike the SupNN (super-ellipse shaped) test 
cases, we have not simulated the mushroom shaped colony for long enough time to reach a quasi¬ 
steady state. Up to time t = 1 s, the mushroom shaped colony exhibits continuous deformation 
characterized by a near-constant IB point displacement velocity of approximately 60 }im / s (not 
shown in the figure). Fig. 13(b) shows the edges in the deformed triangulation after 0.75 seconds, 
colored according to the local strain value. For the sake of clarity, we have only shown the edges 
with strain greater than 0.1 (where positive strain corresponds to a spring that is stretched relative 
to the resting configuration). Clearly, the base of the colony near the substratum experiences the 
highest strain, with values near 0.5. The next highest edge strains (between 0.25 and 0.5) occur 
in sizable portions of the head and base as well as a few points near the midsection of the stem 
region, which can be seen in the zoomed view in Fig. 13(c). An alternate view of the stem is shown 
in Fig. 14(a), with the edges having negative strain colored green, while edges with strain greater 
than 0.25 are colored magenta (and all other edges shown in black). The high strain (magenta) 
edges are clearly aligned along the long axis of the colony, whereas those experiencing compression 
(negative strain, green) are aligned at an angle of 45-60 degrees to the main colony axis. 

We now illustrate a critical drawback of the spring strain-based detachment methodology that 
has not received attention in earlier biofilm IB studies [1,27]. We assume that detachment is ini¬ 
tiated after 0.75 s of deformation and apply a critical strain threshold of e max = 0.25 that is signifi¬ 
cantly lower than the value 1.0 used in these other studies. Performing a single detachment step 
yields the modified spring network shown in Fig. 14(b) (noting that in an actual detachment sce¬ 
nario, these springs would be severed gradually over time instead of all at once). On comparing 
Figs. 14(a) and Fig. 14(b), we see that the edge connectivity around many IB nodes in the stem 
has been altered so that there are now a significant number of rectangular elements in the place 
of triangles. Based on the work of Lloyd et al. [35] we know that two spring networks, one built 
of triangles and the other with rectangles, (but otherwise having the same edge length and spring 
stiffness) will approximate equivalent elastic continua that have different Young's modulus. There¬ 
fore, the spring cutting operation we just performed has effectively introduced an instantaneous 
change in the local mechanical stiffness of the biofilm, which is clearly undesirable and hence is a 
major disadvantage of the spring strain-based detachment strategy. 
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Figure 13: (a) Deformation of the mushroom-shaped biofilm colony at various times in the interval f€[0,l]s. (b) Edges 
colored according to strain at f = 0.75s. (c) zoomed view showing edges inside dotted box in (b). 
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(a) 


Figure 14: (a) Edges in the stem region inside the dotted box i 
0.25 and green indicates negative strain, while all other edges a 
edges with strain greater than 0.25 at 1 = 0.75s. 



(b) 


Fig. 13(b). Magenta identifies edges with strain above 
colored black, (b) Spring connectivity after cutting all 


4.4.2 Equivalent continuum stress-based detachment 

We next apply the equivalent continuum stress-based detachment strategy to the same problem. 
The stress tensor components are computed using Eq. (3.12) after which the von Mises yield stress 
is computed at each IB node using Eq. (3.13). Fig. 15 displays the von Mises stress value inside the 
biofilm region at four time instants during its deformation, where we only color those points with 
stress above the threshold 5 dyne/cm 2 . The von Mises stress has relatively low values throughout 
most of the colony except in three areas: the stem region, near the base where the colony attaches 
to the substratum, and in portions of the biofilm-fluid interface that are subject to large fluid shear. 
The von Mises stress is discontinuous as well as noisy, which is consistent with other numerical 
simulations of microstructural stress inside granular materials [6,20]. However, we emphasize that 
this behaviour is in stark contrast with the apparently smooth von Mises stress field obtained by 
Towler et al. [54] for a simpler biofilm shape using a finite element simulation. 

We make no attempt here to draw any explicit correspondence between results from the strain- 
and stress-based detachment strategies because the spring strain-based threshold parameter e max 
cannot be translated into a von Mises yield stress value. However, we can still compare the two by 
simulating detachment for the equivalent continuum stress methodology using the same param¬ 
eters. This time, we initiate detachment after 0.25 s of deformation and choose the two distance 
threshold parameters e suh = e ext = 2 urn (refer to Section 3.3). This divides the mushroom shaped 
colony into three zones: in zone 1 near the substratum, the adhesive strength is E fl d/j=10 dyne/cm 2 ; 
for zone 2 near the biofilm-fluid interface, the interfacial cohesive strength is = 0.1 dyne/cm 2 ; 

and for zone 3 in the interior, the cohesive strength is = 2.5 dyne/cm 2 . Based on these pa¬ 
rameter values. Fig. 16 depicts those IB points that will be severed according to the local value 
of von Mises stress, with the red points indicating detachment in zone 2, while the green points 
corresponding to detachment in zones 1 and 3. From this figure it is clear that two regions within 
the colony neck experience complete rupture, which is in agreement with what one would expect 
from solid mechanics. Furthermore, the detachment criterion is active along the entire biofilm- 
fluid interface, which corresponds to an erosion process. 

In summary, our new detachment strategy based on equivalent continuum stress provides an 
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Figure 15: Von Mises stress inside the deformed biofilm at a sequence of times. For purpose of clarity, only points where 
the stress exceeds 5dyne/crrr are colored. 
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Figure 16: Points where detachment occurs in the biofilm colony at time t= 0.25s (zoomed-in view on the right). Detachment 
parameters are e sub = e ext = 2, Z adh = 10, E'^-2.5 and E^ = °.!. 


unambiguous method for performing biofilm detachment that is also consistent with methods 
employed by other continuum biofilm models. The detachment of any IB point proceeds by cutting 
of all springs attached to it so that the topology of the spring network remains triangular, thereby 
avoiding a major disadvantage of the spring strain-based detachment strategy However, this 
approach does introduce some additional computational work in each time step for evaluating the 
equivalent continuum stress and distance functions at each IB node, not to mention maintaining 
all of the relevant data structures. 


4.5 Biofilm deformation and internal stress 

This section presents a final series of simulations that study the biofilm stress distribution in more 
detail and draw specific conclusions regarding the various modes of detachment (sloughing or 
erosion). We also provide evidence that questions the validity of another class of detachment 
models based on a detachment speed function. 

The steady-state von Mises stress for the SEMI20 and SUP75 simulations from the previous 
section are depicted in Fig. 17 for two values of biofilm spacing, = 50 and 400 pm. In all cases, 
the stress is lowest in the interior of the biofilm and largest along the biofilm-fluid and biofilm-wall 
interfaces, with the absolute maxima occurring on the wall near the leading and trailing corners. 
We also observe that the proportion of the biofilm experiencing high von Mises stresses increases as 
the spacing parameter Dj, increases, which is consistent with the results from Section 4.3. Finally, 
the stresses experienced in a long, thin colony such as SUP75 are significantly larger (note the 
increase in the colormap scale by a factor of 10 for plots in the bottom row in Fig. 17). 

The corresponding plots of von Mises stress along the bottom edge of the colony are shown 
in Figs. 18(a,b) from which it is clear that the stress peaks at the up/downstream corners. These 
plots actually terminate a short distance away (5-8 IB points) from the corners because the stress 
maxima at the corners is 8 to 10 times the value in the interior. We do the same for the stress along 
the biofilm-fluid interface in Figs. 18(c,d), from which we see that the interfacial shear stress attains 
a local maximum near the biofilm tip where the fluid shearing force is largest. On comparing the 
stress plots for the wall and biofilm-fluid interfaces, it is clear that increasing the colony spacing 
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Figure 17: Von Mises yield stress at steady state for cases SEMI20 and SUP75, with K b j o =0.75 and D b =5 0 and 400//m. For 
ease of visualization, only IB points with stress above a threshold of 0.050 dyne/cm 2 in (a) and (b) and 0.20 dyne/cm 2 
in (c) and (d) are shown. 
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(b) SUP75, biofilm-wall interface 




(d) SUP75, biofilm-fluid interface 


Figure 18: Von Mises yield stress adjacent to the wall (a,b) and on the biofilm-fluid interface (c,d) for the same simulations 
as in Fig. 17. 
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D/, has a much more pronounced effect on the biofilm-fluid interface than on the wall stresses. We 
conclude from these results that for widely-spaced colonies, detachment is more likely to occur by 
surface erosion than by sloughing from the substratum. This behavior is to be expected since we 
have used a relatively large value of the spring stiffness connecting the wall and biofilm (compared 
to the value inside the biofilm). However, if we were to use a lower value to mimic weak wall 
attachment then we could capture the competition between surface erosion and sloughing modes 
of detachment, with the dominant mode being determined by the relative sizes of the adhesive 
(E adf,) and exterior cohesive (E'’^) stress thresholds. 

Upon closer investigation of the plots in Fig. 17(c,d), we note that for case SEMI20 the von Mises 
stress increases along the entire exposed biofilm surface as D\, increases, whereas in case SUP75 
only the left face experiences an increase. This suggests that surface erosion is more uniform for 
circular colonies, whereas long thin colonies will tend to erode only along the upstream face. Based 
on this result, and the fact that von Mises stresses are much larger for SUP75, it seems reasonable 
to suppose that the enhanced surface erosion observed in elongated colonies could be a precursor 
to formation of streamers that are observed in experiments at low Reynolds number [46,56]. An 
in-depth investigation of streamer formation is beyond the scope of this work but could form the 
basis for an interesting future study based on our IB model. 

Additional insight can be gained by comparing plots of interfacial shear stress computed earlier 
in Figs. 12(a,d) (with k f);o = 0.75) with the corresponding von Mises stress curves for the same 
biofilm colonies in Figs. 18(c,d). For the SEMI20 case, the interfacial fluid stress and von Mises 
stress curves have the same general shape along the central portion, with the main difference 
being that the von Mises stress increases towards the corners whereas the interfacial shear stress 
does not. For the elongated colony in SUP75, both stresses are asymmetric about the colony apex 
(near IB point 360) and although both experience a rapid decrease to the right (downstream) of the 
apex, the upstream behaviour is very different. In particular, the interfacial fluid stress increases 
gradually on the left toward the apex, while the von Mises stress sustains a relatively large value 
on the left with a more rapid rise to the maximum. 

These differences just mentioned point to an important error in another commonly-used de¬ 
tachment approximation based on a detachment speed function. In this approach, rather than ex¬ 
plicitly solving the continuum equations for mechanically-induced detachment, they account for 
these effects instead by specifying a local speed at which the biofilm-fluid interface recedes into the 
biofilm [36,65]. For example, this detachment speed function may depend on the local interfacial 
shear stress [13] or interfacial curvature [65]. Based on the IB results above in which significant 
differences occur between fluid shear stress and von Mises stress along the biofilm-fluid interface, 
it is clear that even for a reasonably stiff biofilm a detachment speed function depending on in¬ 
terfacial shear stress is incapable of correctly capturing detachment dynamics at all points along 
the interface. It is possible that an alternate speed function could be found that accounts for the 
variation in von Mises stress along the interface that has been identified in this study and so this 
would be an worthwhile subject for further study. 

As a further illustration of the detachment process. Fig. 19 shows the portion of the biofilm that 
will detach at steady state for cases Sup25, Sup50 and Sup75 with kj, !O =0.75 and Dj,= 50 or 400f<m. 
The remaining parameters in the detachment algorithm from Section 3.3 are e suh = e ext = 2 pm, 
E e ff h = 0.1 dyne/cm 2 , Ydff = 1 dyne/cm 2 and E (;rf /, = 5 dyne/cm 2 . In practice, the removal of IB 
points by detachment will alter the forces acting on the colony, which in turn induces further de¬ 
formation; this process repeats until no further detachment is possible. Simulations implementing 
this alternating detachment/deformation process were conducted for a few selected cases; how¬ 
ever, the time step limitation in our IB algorithm (At « 10 5 s) precluded integrating the solution 
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(a) SUP25, D(, = 50/rm 




(e) Sup 75, Df, = 50/;m 





Figure 19: Plots indicating the portion of the biofilm colony (at steady state) that will detach for cases SUP25, SUP50 
and SUP75 with Kj, io = 0.75. Sub-figures (a,c,e) have colony spacing Dj, = 50/(m, while (b,d,f) have Dj,=400/rm. Other 
parameters: e sub = e ext = 2}irt\, = 0.1 dyne/cm 2 , = l dyne/cm 2 and Z„ rf ^ = 5dyne/cm 2 . 


over the long time intervals required. Consequently, the scope of this work is restricted to intro¬ 
ducing the equivalent continuum stress based detachment strategy and validating the results on a 
range of colony shapes. Efforts to develop a more efficient implementation for a complete defor¬ 
mation and detachment strategy are currently underway with a more efficient IB algorithm and 
will form the basis for a future publication. 


5 Conclusions 

We employed a 2D immersed boundary method to simulate the deformation of a periodic array of 
uniformly-spaced, wall-bounded, weak biofilm colonies in response to a linear shear flow. In order 
to capture different stages of biofilm growth under mass transfer-limited conditions, we chose a 
family of biofilms having the same generic shape (sections of a super-ellipse) but with increasing 
aspect ratios (fixed width, increasing height). Actual biofilm colonies behave mechanically like vis¬ 
coelastic solids and we mimic this behaviour by replacing the biofilm with a network of Hookean 
springs corresponding to the edges in a quasi-uniform triangulation of the colony. 

We began by performing a parametric study that investigated the effect of colony spacing and 
spring stiffness on the drag/lift forces and interfacial shear stress acting on the biofilms. The main 
results of this parametric study can be summarized as follows: 

• Varying spring constant: At low shear rates, colonies with even a moderate spring stiffness 
of Kjjio ^ 10 are able of resisting large deformation forces. Weak biofilm colonies with k/, ;o ^ 1 
experience increased drag and larger deformations with a maximum displacement in the 10's 
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of microns. These larger deformations are accompanied by a change in the interfacial shear 
stress profile wherein stress increases along the upstream face and decreases downstream. 

• Varying colony spacing: For low to moderate biofilm stiffness (k/, ;o < 10) reducing the colony 
spacing from 400 to 50 pm reduces the drag by as much as 50-100%, accompanied by a 
change in the shear stress profile along the biofilm-fluid interface. The interfacial shear stress 
in colonies with large aspect ratio (SUP75) differs significantly from ones with small aspect 
ratio (SEMI20). 

• Varying colony shape: It is possible for biofilm colonies to grow into tall structures with large 
aspect ratio even if they are weak mechanically (ic^o < 10) owing to the protection from sur¬ 
face erosion afforded by being in close spatial proximity to other colonies. We believe that 
this result will carry through to 3D (although to a lesser degree) and this is an issue that we 
plan to investigate further in a future study. 

We also developed a new method for initiating biofilm detachment using averaged equivalent 
continuum stress, which we implemented within our IB framework. We overcame problems en¬ 
countered in other spring strain based detachment strategies such as in [27] and [1], The following 
conclusions can be drawn regarding detachment: 

• Variation in von Mises stress: Increasing the spacing between colonies leads to an increased 
tendency for surface erosion instead of sloughing from the wall. Based on the von Mises 
stress along the biofilm-fluid interface, we concluded that semi-circular biofilm colonies will 
undergo roughly uniform surface erosion, while colonies with larger aspect ratios erode pre¬ 
dominantly along the upstream face. 

• Correlation between fluid stress and von Mises stress: The fluid stress and von Mises yield stress 
along the biofilm-fluid interface differ substantially. We conclude that biofilm dynamics 
based on a detachment speed function approach [13,36] (where detachment speed is a func¬ 
tion only of local interfacial fluid stress) cannot capture the actual detachment behaviour re¬ 
sulting from excessive straining. This highlights the importance of using detachment strate¬ 
gies that accurately capture the biofilm mechanics. 

The main advantage of our detachment algorithm is that it provides a uniform framework 
for handling biofilm deformation, surface erosion and sloughing through the use of a continuum 
mechanics-based detachment model that employs measured biofilm mechanical properties. The 
primary disadvantage of our approach as implemented herein is the high computational cost; in 
particular, the small time step required for stability reasons combines with the extra work of simu¬ 
lating detachment to make long-time computations of simultaneous deformation and detachment 
impractical. We emphasize that this is not a limitation of the immersed boundary approach, but 
rather our specific implementation that uses a simple explicit time-stepping strategy. We conclude 
therefore that the full potential of our detachment algorithm can only be realized in combination 
with either a (semi-)implicit time-stepping approach or an efficient parallel implementation (such 
as [62]), which will be the subject of future work. 
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